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1.  INTRODUCTION 


Because  of  their  relevance  to  the  aerodynamics  of  vertical-takeof f-and- 
landing  aircraft  in  ground  effect,  considerable  research  has  recently  been 
devoted  to  isolated  and  interacting  impinging  jets.  Figure  1  illustrates 
these  configurations  schematically  in  planar  view.  For  the  case  of  a  single 
jet.  Figure  1(a),  there  are  four  simultaneous  flow  regimes  of  differing 
character:  the  transition  region  in  which  the  jet  velocity  profile  changes 

from  a  uniform  to  a  fully  developed  distribution;  the  fully  developed  region; 
the  impingement  region  in  which  the  fluid  is  deflected  by  the  ground  plane; 
and  the  wall-jet  region  in  which  the  flow  spreads  radially  from  the  stagnation 
point.  For  the  case  of  two  interacting  jets,  Figure  1(b),  these  four  regimes 
exist  for  each  of  the  impinging  jets,  and,  in  addition,  a  fountain  region  is 
formed  through  the  collision  of  the  opposing  wall  jets. 

Although  the  flow  patterns  illustrated  in  Figure  1  appear  simple,  they 
are  nevertheless  difficult  to  compute  because  they  are  characterized  by 
significant  turbulence  levels,  strong  pressure  gradients,  and  an  elliptic 
character  which  precludes  analysis  by  parabolic  marching  schemes.  The  ground- 
plane  pressure  distribution  in  the  impingement  zone,  which  is  governed  by  the 
dynamics  of  the  flow,  can  be  predicted  adequately  by  inviscid  analytical 
methods,  but  calculation  of  entrainment  effects  and  velocity  profiles  in  the 
wall  jet  and  fountain  necessitates  the  inclusion  of  turbulence. 

For  isolated  impinging  jets,  both  two-  and  three-dimensional,  Rubel  has 
solved  the  conservation  equations  for  inviscid  flow  in  finite-difference  form, 
References  1  and  2.  These  analyses  provided  good  agreement  between  measured 
and  computed  ground-plane  pressure  distributions  for  both  normal  and  oblique 
impingement.  For  two-dimensional  (planar)  impinging  jets,  Agarwal  and  Bower, 
Reference  3,  included  viscous  effects  by  solving  the  time-averaged  Navier- 
Stok.es  equations  for  both  incompressible  and  compressible  flow  and  using  the 
two-equation  Jones-Launder  model  to  obtain  closure.  This  work  yielded  good 
agreement  with  pressure  and  velocity  data  and  provided  information  on  the 
variation  of  turbulence  length  scale  with  Reynolds  number. 

For  two  interacting,  Impinging  jets  with  fountain  formation,  most  tech¬ 
niques  developed  to  predict  these  flows  use  the  so-called  "modular  analyses," 
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such  as  those  described  by  Kotansky  et  al.,  Reference  4,  and  Siclari  et  al.,  j 

Reference  5.  With  this  approach,  the  various  regions  of  the  flowfield,  such 
as  the  free  jets  and  the  wall  jet,  are  modeled  using  a  separate  system  of  \ 

equations  for  each  zone,  which  are  generally  based  on  integral  conservation- 
law  constraints.  The  computed  properties  are  then  iteratively  pieced  together  i 

to  provide  the  quantities  of  interest.  This  method  works  well  for  relatively 
wide  jet  spacing  (S/D  >  6)  where  the  fountain  properties  are  dependent  pri¬ 
marily  on  the  properties  of  the  approaching  wall  jets.  However,  as  S/D  is 
reduced  below  6,  there  is  a  strong  interaction  between  the  impinging  jets  and 
the  fountain,  as  established  in  the  flow  visualization  work  of  Sarlpalli, 

Reference  6.  The  coupling  between  the  primary  jets  and  the  fountain  invali¬ 
dates  any  prediction  model  that  computes  the  upwash  based  only  on  the  wall-jet 
properties.  Siclari  et  al.,  Reference  7,  have  formulated  a  "modular  analysis" 
for  closely  spaced  jets,  but  this  scheme  is  restricted  by  the  need  for  numer¬ 
ous  empirical  constants  which  are  configuration-dependent. 

A  more  general  approach  for  computing  dual-jet  impingement  flows  is  the 
solution  of  the  time-averaged  Navler-Stokes  equations  in  three  dimensions. 

This  problem  is  considerably  more  complex  than  the  planar  impinging-jet 
problem  for  the  following  two  reasons: 

(1)  Computer  storage  limitations  restrict  the  number  of  finite-difference 
grid  points  available,  which  limits  the  maximum  Reynolds  number  for 
which  convergent  or  accurate  solutions  can  be  obtained. 

(2)  Even  with  coarse  grids,  the  computer  time  necessary  to  achieve  a 
three-dimensional  solution  is  at  least  an  order  of  magnitude  greater 
than  the  time  necessary  to  compute  a  two-dimensional  solution,  which 
translates  directly  into  cost  considerations. 

Initial  work  at  MDRL  on  solution  of  the  three-dimensional,  time-averaged, 
Navier-Stokes  equations  was  based  on  writing  the  conservation  equations  in 
terms  of  scalar  and  vector  potential  functions  and  vorticity,  using  the  one- 
equation  Glushko  model  to  represent  the  turbulence,  applying  experimentally 
determined  boundary  conditions  in  the  near  field,  and  solving  the  transport 
equations  using  Hoffman's  augmented-central-dif ference  algorithm  (Reference 
8).  With  this  approach,  flowfields  were  computed  for  equal-  and  unequal- 
strength  jets  with  normal  and  oblique  impingement,  and  representative  results 
given  in  the  form  of  particle  pathline  traces  are  illustrated  in  Figure  2. 
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Figure  2.  Streamline  plots  for  equal-  and  unequal-strength  jets  with  normal  impingement 
(Re  =  VjcejD2/M  *  100). 


Shortcomings  of  this  formulation  are  the  need  to  specify  in  the  Glushko 
model  the  length-scale  distribution  of  the  turbulence,  which  is  a  complex 
variation  for  jet-impingement  flows,  and  the  need  to  establish  the  boundary 
conditions  in  the  near  field,  which  generally  must  be  evaluated  from  experi¬ 
mental  results.  To  overcome  these  limitations,  the  following  are  the  specific 
objectives  of  the  current  work  on  solution  of  the  time-averaged  Navier-Stokes 
equations  for  three-dimensional  turbulent  impinging- jet  flows: 


(1)  Utilize  the  Jones-Launder ,  two-equation  turbulence  model,  thereby 
eliminating  the  need  to  specify  the  length-scale  distribution. 

(2)  Utilize  coordinate  transformations  which  permit  imposing  the  boundary 
conditions  in  the  far  field  where  relatively  simple  constraints  can  be 
placed  on  the  flow  variables,  thereby  eliminating  the  need  to  intro¬ 
duce  empirically  based  boundary  conditions  in  the  near  field. 

(3)  Solve  the  governing  equations  directly  for  the  velocity  and  vorticity, 
thereby  eliminating  the  need  for  scalar  and  vector  potential  functions 
which  require  additional  computer  storage. 

(4)  Discretize  the  transport-type  equations  using  a  third-order-accurate 
upwind-difference  scheme  which  is  more  efficient  than  the  augmented 
central-difference  algorithm. 

This  report  describes  the  governing  equations,  the  numerical  solution 
technique,  and  computed  jet-impingemc rt  flowfields.  The  latter  include  iso¬ 
lated  jets  with  both  normal  and  oblique  impingement  and  interacting  equal-  and 
unequal-strength  jets  with  fountain  formation. 


2.  THE  FLOWFIELD  MODEL 


This  section  presents  the  governing  equations  used  to  describe  the  turbu¬ 
lent,  incompressible,  steady  flow  of  the  three-dimensional  impinging-jet  con¬ 
figurations.  The  complete  time-averaged  Navier-Stokes  equations  and  Jones- 
Launder  turbulence  model  equations  are  given  first  in  terms  of  the  velocity 
components  and  pressure,  and  then  they  are  rewritten  in  terms  of  velocity  and 
vorticity.  Coordinate  transformations  are  applied  to  the  equations,  and  the 
types  of  boundary  conditions  imposed  in  the  analysis  are  defined. 

2.1  Governing  Equations 

Considering  a  space-fixed  reference  through  which  the  fluid  moves,  the 
time-averaged  continuity  and  momentum  (Navier-Stokes)  equations  for  steady 
incompressible  flow  are  given  below  in  tensor  notation. 


Conservation  of  mass: 


la 

3x. 


0 


(1) 


Conservation  of  momentum  in  the  ic^  direction: 
3(u  .u, ) 

+  ^ 


— j- -1—  +  iE - L.  (f  _  pu  'u  ') 

3Xj  3x1  3Xj  VTij  pui  j  ' 


(2) 


The  conventional  notation  is  used  in  which  p  is  the  density,  u^  the  velocity 

f  h 

component  in  the  iu  direction,  p  the  static  pressure,  and  x .  .  the  shear 
stress  in  the  i*"  direction  on  a  surface  normal  to  the  jcn  direction. 


The  previous  equations  were  derived  using  Reynolds  decomposition  in  which 
a  general  flow  variable  $  is  expressed  as  the  sum  of  a  time-averaged  compo¬ 
nent  <t>  and  a  fluctuating  component  $',$=$+$'.  The  result  is  a  system  of 
time-averaged  conservation  equations  which  have  the  same  form  as  the  instan¬ 
taneous  conservation  equations  with  the  exception  of  the  Reynolds  stress 
terra  -  pu^ 'u  ^ 

In  order  to  obtain  a  closed  system  of  equations,  the  Reynolds  stress  term 
must  be  related  to  the  time-averaged  quantities  that  describe  the  mean  flow. 
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The  present  work  uses  the  high-Reynolds-number  form  of  the  two-equation  Jones- 
Launder  (k-e)  turbulence  model.  Reference  9.  With  this  approach,  the  effec¬ 
tive  shear-stress  term  in  Equation  (2)  is  represented  by 


_  i 

Tlj  '  pui'uj'  '  “eff  l nj  +  3^1  '  f  5ljpk>  <3> 

where  Is  the  effective  viscosity  (the  sum  of  the  molecular  and  turbulent 

components,  ■  y  +  ut)  and  k  is  the  turbulent  kinetic  energy.  In  the 

high-Reynolds-number  form  of  the  Jones-Launder  turbulence  model,  is  related 
to  the  turbulent  kinetic  energy  k  and  its  kinematic  rate  of  dissipation  c  by 


where  c  is  a  constant  with  a  value  of  0.09.  In  this  turbulence  model,  k 
U 

and  e  are  computed  from  the  transport  equations  given  below  in  tensor 
notation. 


Transport  equation  for  turbulent  kinetic  energy: 


-  3k 

PUj3^ 


JL  L  5!s_\  +  u  +  !!i. 

j  3xi  \k3xi)  c  y*i  5xi /  3xj 


Transport  equation  for  turbulent  dissipation: 

In  Equations  (5)  and  (6)  the  diffusion  coefficients  y,  and  y  are  defined  by 
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The  constants  appearing  in  the  turbulence-model  transport  equations  have  the 
following  values:  =  1.44,  C2  =  1.92,  X^  =  1.0,  and  X£  =  0.9091. 

Although  the  governing  equations  are  written  in  velocity-pressure  form, 
they  are  solved  in  velocity-vorticity  form  to  simplify  the  solution  procedure 
by  olirainating  the  pressure  gradient  from  the  momentum  equations.  The 
vorticity  is  defined  in  tensor  notation  by 


$2.  -  e.  ..  ,  (9) 

1  ijk  3x^ 

where  e.,,  is  the  alternating  tensor.  Taking  the  curl  of  Equation  (2)  and 
ijk 

utilizing  Equation  (9)  results  in  a  transport  equation  for  the  vorticity.  A 
Poisson  equation  for  the  velocity  field  is  derived  by  combining  the  relation 
for  the  conservation  of  mass,  Equation  (1),  with  Equation  (9). 

Completing  the  steps  outlined  in  the  previous  paragraph,  introducing  a 

characteristic  length  D  and  characteristic  velocity  VQ  to  normalize  the 

equations,  and  expanding  the  equations  in  cartesian  coordinates  (x,y,z)  with 

velocity  components  (u,v,w)  and  vorticity  components  (fl  ,JJ  )  results  in  the 

x  y  z 

system  of  dimensionless  equations  given  below. 


Poisson  equation  for  u: 


Poisson  equation  for  v: 


Poisson  equation  for  w: 


Transport  equation  for 


/a2*}  a2: 
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\3x  3y 
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Transport  equation  for  Qy: 


(1  +  Re 
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Transport  equation  for  e: 


(17) 


The  only  dimensionless  parameter  that  enters  the  system  of  equations  is 
the  Reynolds  number,  Re  -  VqD/v.  The  location  in  the  flowfield  where  VQ  and  D 
are  chosen  will  be  defined  subsequently  with  reference  to  the  specific 
impinging- jet  configurations. 


2.2  Coordinate  Transformations 

In  the  present  analysis,  a  transformation  is  used  in  each  of  the  three 
coordinate  directions.  This  approach  permits  mapping  a  relatively  restricted 
computational  domain  into  a  larger  physical  domain  and  distributing  points  in 
the  finite-difference  solution  of  the  equations  in  regions  of  the  flowfield 
where  property  gradients  are  the  largest. 

The  transformed  coordinates  are  denoted  by  an  overbar:  x  *  x(x), 
y  “  y(y),  and  z  ■  z(z).  It  follows  that  the  required  derivatives  of  a  general 
flow  variable  4  with  respect  to  the  physical  coordinates  are  related  to  the 
derivatives  with  respect  to  the  transformed  coordinates  by  the  expressions 
given  below. 


H 


(18) 


dx  3$ 

3x  *  dx  - 
3x 

(18) 

24-  (§) 

2  2  2- 
|  3 ^  d  x  34> 

1  _o  '  2  — 

(19) 

a*2  V*) 

3x  dx  3x 

3<t>  dy  3<t> 

ay  '  3j 

(20) 

aii  .  I ii'1 

ay2  H 

2  2  2- 
|  +  ,d— y  14. 

1  -2  2- 
3y  dy  3y 

(21) 

34  dz  3* 

(22) 

d24  /dz ' 

2  “  Idz 
3z  \  / 

i  2  2  2- 

l  3  4  +  d  z  34 

1  3 z2  dz2  3z 

(23) 

The  Poisson  equations  for  the  velocity  components.  Equations  (10)  through 
(12),  can  be  written  in  the  following  general  form  in  terms  of  the  physical 
coordinates : 


ifi  +  i±  +  »i| .  „ 

3x  3y  3z 


♦  ’ 


(24) 


where  $  -  u,  v,  or  w  and  a ^  denotes  the  corresponding  source  term.  In  terms 
of  the  transformed  coordinates.  Equation  (24)  assumes  the  following  form: 


1  3xZ  2  3x  1  3y2  2  3y  3* 


3*  - 

-  -*■  -  a.  , 

2  3z  ♦ 


(25) 


where 


dx 

dx 


d2; 


>  “2  ■  ZJ  •  81  '  §  ’  '2  '  ^2  •  T‘  '  &  ’  and  ‘  ^2  • 


dx 


The  overbar  on  the  source  term  denotes  that  the  derivatives  it  contains 
have  been  rewritten  in  terms  of  the  transformed  coordinates  according  to 
Equations  (18)  through  (23). 


Similarly,  the  transport  equations  for  the  vorticity  components  and  the 
turbulence  quantities.  Equations  (13)  through  (17),  can  be  written  in  terms  of 
the  physical  coordinates  in  the  following  general  form: 


where  $  *  fl  ,  8  ,  n  ,  k,  or  e.  The  corresponding  coefficients  in  the  respec- 
x  y  z 

tive  equations  are  denoted  by  a,,  8.,  y.,  and  6  ,  and  the  source  term  is 

?  v  ♦  v 

denoted  by  a  .  When  Equation  (26)  is  written  In  terras  of  the  transformed 

♦ 

coordinates,  the  following  equation  results: 


V?  §  +  V*  +  V?  $  +  (v*  +  Vi) 


M 

3x 


+  (v**Vi)  fy+  (V2  +  Vi) 


■♦r**Vi)  it'5*  • 


(27) 


where  Oj,  a 2#  8^,  $2»  Yj»  and  y2  maintain  the  same  definitions  introduced  with 
reference  to  Equation  (25). 

The  analytic  functions  used  to  define  the  coordinate  transformation 
derivatives  are  different  for  the  various  impinging-jet  configurations  con¬ 
sidered  and  are  defined  in  a  subsequent  section. 

2.3  Boundary  Conditions 

In  the  system  of  equations  presented  in  Section  2.1,  the  eight  elliptic 

partial-differential  equations  must  be  solved  for  the  following  eight  flow 

properties:  u,  v,  w,  ft  ,  ft  ,  (2  ,  k,  and  e.  The  manner  in  which  the  boundary 

x  y  z 

conditions  are  Imposed  on  these  equations  depends  on  the  nature  of  the 
boundary  plane  in  the  physical  domain  where  the  impinging-jet  flowfleld  is 
solved. 

Consider  first  the  inflow  plane  where  one  or  two  jets  enter  the  solution 
domain  a  distance  y  ■  K  above  the  ground.  The  velocity  and  turbulence 
property  profiles  throughout  the  entering  jets  are  assumed  to  be  known: 
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u  ■  u(x,  h,  z),  v  -  v(x,  h,  z),  w  -  w(x,  h,  z),  k  -  k(x,  h,  z), 
and  e  ■  e(x,  h,  z). 


(28) 


The  profiles  of  the  vorticity  components  through  the  entering  jet  follow  from 
the  specified  velocity  distributions  and  the  defining  equations  for  the 
vorticity  components. 


a 


X 


3w  3v 

3y  3z 


(29) 


n 


y 


3u 

3z 


3w 

3x 


(30) 


n 


z 


3v  3u 
3x  3y  * 


(31) 


On  a  solid  surface  (either  the  impingement  plane  or  the  blocking  plate 
through  which  the  jet  discharges),  the  velocity  and  turbulent-kinetic-energy 
boundary  values  follow  from  the  no-slip,  impermeable-wall  condition. 


u=*v”w  =  k  =  0. 


(32) 


The  turbulent  dissipation  is  specified  by  taking  the  limiting  form  of  Equation 
(16)  at  the  wall. 


e 


1 

Re 


(33) 


where  n  denotes  the  coordinate  normal  to  the  surface.  The  vorticity  profiles 
follow  from  Equations  (29)  through  (31). 

On  an  outflow  plane  where  the  wall  jet  exits  the  solution  domain,  it  is 
assumed  that  there  are  no  gradients  with  respect  to  the  coordinate  normal  to 
the  outflow  surface.  (The  outflow  planes  are  taken  sufficiently  far  from  the 
impinging  jet  centerline  such  that  this  assumption  is  valid.) 


3u  _  3v  _  3w  _  _  5.  .  5i  .3Ji.il 

3n  3n  3n  "  3n  3n  3n  3n  3n 


(34) 


Finally,  the  remaining  boundary  surface  to  consider  is  a  symmetry  plane 
passing  through  the  impinging  jet  or  fountain.  On  this  plane  the  normal 
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1 


I 


I 

I 


velocity  component  vanishes,  and  the  derivatives  of  the  other  velocity 
components  with  respect  to  the  normal  coordinate  are  zero. 


u 

n 


*  0  . 


I  l#n 

The  normal  gradients  of  the  turbulence  quantities  vanish  on  the  symmetry 
plane ,  * 


(35) 


3k  _  3e 
3n  ”  3n 


0  , 


(36) 


and  the  vorticlty  values  follow  from  their  defining  equations. 


Specific  solution  domains  for 
responding  boundary  conditions  are 


the  lmplnging-jet  geometries  and  the  cor- 
presented  in  Section  4. 


» 


15 


3.  THE  NUMERICAL  SOLUTION  SCHEME 


This  section  is  a  discussion  of  the  finite-difference  technique  used  to 
solve  the  partial-differential  equations  that  describe  the  flowfield.  The 
Poisson-type  equations  are  discretized  using  the  standard  central-difference 
algorithm,  and  the  transport-type  equations  are  discretized  using  a  third- 
order-accurate  upwind-difference  scheme.  The  discretization  of  the  boundary 
conditions  and  the  iterative  solution  of  the  coupled  system  of  equations  are 
described. 

3.1  Discretization  of  the  Poisson-Type  Equations 

As  given  in  Section  2,  the  Poisson  equations  for  the  velocity  components 
can  be  written  in  the  form 


2  3  <j>  .  36  .  „2  3  6  .  36  .  2  3  6,  36 

,  — L  +  a,  -*■  +  g.  —f  +  0,  +  y .  — £  +  Y, 

3xz  z  3x  3y  L  By  3z  3z 


(37) 


where  6  ■  u,  v,  or  w. 

To  write  Equation  (37)  in  finite-difference  form,  the  three-dimensional 
nodal  network  shown  in  Figure  3  is  Introduced.  Since  coordinate  transforma¬ 
tions  are  applied  to  the  differential  equations,  a  nonuniform  finite- 
difference  grid  is  not  required  to  distribute  nodes  in  regions  of  large 


I  =  maximum  i  index 


Figure  3.  The  three-dimensional  finite-difference  stencil. 


s 

f 

i 


property  gradients.  Consequently,  a  uniform  grid  spacing  h  is  used  in  each  of 
the  transformed  coordinate  directions  so  that  x  =  h(i  -  2),  y  *  h(j  -  2),  and 
3  =  h(k  -  2),  where  i,  j,  and  k  are  the  respective  node  indices.  The  bound¬ 
aries  of  the  computational  domain  are  the  planes  i  -  2,  j  *  2,  k  =  2,  i  ■  I, 
j  =  J,  and  k  -  K;  the  planes  i  =  1,  j  =  1,  k  *  1,  i  -  1+  1,  j  »  J  +  1,  and 
k  *  K  +  1  are  exterior  planes  required  in  discretization  of  the  transport-type 
equations. 

The  conventional  central-difference  algorithm  is  introduced  with  the 
derivatives  centered  at  the  interior  point  (i,j,k).  The  approximations  for 
the  first  derivatives  are 


3$ 

3x 


i+l,j,k  ~  *i-l.j.k 
2h 


i.j.k 


a<p 

ay 


♦l,J-H,k  " .♦lt.J-l,k 
2h 


i.j.k 


it 

3z 


♦i,J,k-fl  “  ♦i.J.k-l 
2h 


i.j.k 


and  the  approximations  for  the  second  derivatives  are 


08) 


(39) 


(40) 


ifi 

a;2 


i!± 

ay2 


.2 


.  ♦i,j+i,k  ~  2*i,.i,k  +  ♦i.j-i.k 

.2 


(41) 


(42) 


i.j.k 


ift 

3z2 


-  ♦l.j.k-t-l  ~  2»i,j,k  +  »i,j,k-l  >  (43) 

i,j,k  h2 

When  Equations  (38)  through  (43)  are  substituted  into  Equation  (37),  the 
discretized  form  of  the  Poisson  equation  is 
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Vi+l,j,k  +  Vi-1,  .i,k  +  Vi,j+l,k  +  Vi,j-l,k  +  Vi,j,k+1 


+  Vi,j,k-1  +  ¥i,j,k  =  V 


(44) 


The  coefficients  are  defined  by  the  following  relations: 

=  a\  +  0.5ha2 
Pj  *  -  0.5ha2 

P4  -  +  0.5hP2 

p5  ■  B j  -  O.5h02 

P6  ”  Y1  +  °*5h^2 

P7  ”  Y1  _  °*5h^2 
p2  *  ~  2(aJ  +  +  Yj) 

P8  ■  h%- 

To  solve  the  system  of  algebraic  finite-difference  equations, 
(44)  is  rewritten  for  3<j<J-lin  the  following  form: 


(45) 

(46) 

(47) 

(48) 

(49) 

(50) 

(51) 

(52) 

Equation 


Vi,J+l,k  +  Vi,j,k  +  Vi,j-l,k 


P8  "  Vi+l,j,k 


Vi-1, j,k  ‘  Vi,j,k+1  ‘  Vi,j,k-1  * 


(53) 


Applying  Equation  (44)  for  j  ■  3,  $  _  is  known  from  the  boundary  value  and 

i  f  Z  |  K 

is  transferred  to  the  right  side  of  the  equation.  Therefore, 


Vi,4,k  +  Vi,3,k  “  P8  "  Vi+l,3,k  “  Vi-1, 3, k  “  Vi,2,k 
"  Vi,3,k+1  "  *Vi,3,k-l  ’ 


(54) 
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Similarly,  when  Equation  (44)  is  written  for  j  =  J  -  1,  fk^8  known  ^ro® 
the  boundary  value  and  is  transferred  to  the  right  side  of  the  equation. 


P2*i,J-l,k  +  P5*i,J-2,k  P8  Pl*i+l,J-l,k  “  *3*1-1 , J-l ,k 

P4*i, J,k  “  P6*i,J-l,k+l  "  P7*i,J-l ,k-l  *  (5' 

Equations  (53)  through  (55)  are  solved  using  the  tridiagonal  algorithm. 
Reference  10.  Equation  (53)  is  rewritten  in  the  form 

p5d»j-i +  Vj  +  Vj+i  ■  \  *  <51 

where  the  subscript  has  been  added  to  the  P  coefficients  to  indicate  the  j 
node  at  which  they  are  evaluated;  the  i  and  k  subscripts  have  been  dropped 
from  the  unknown  4  values  for  simplicity;  and  Rp^  denotes  the  right  side  of 
Equation  (53).  A  recursion  relation  is  assumed. 


♦j  ‘  qJ  "  ’ 

and  substitution  of  Equation  (57)  into  Equation  (56)  yields 


Y,  “  P5jqj-1 


j  P2j  "  P5jbj-1  P2j  "  P5jbj-1 


From  a  comparison  of  Equations  (57)  and  (58),  it  follows  that 


Yj  "  P5jqj-1 
lj  P2j  "  P5jbj-1 


J  P2j  "  P5jbj-1 
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To  initialize  the  solution,  Equation  (56)  is  written  for  j  =  2, 


(61) 


According  to  the  recursion  formula, 
*2  “  q2  '  b2*3  * 


and  it  follows  from  a  comparison  of  Equations  (61)  and  (62)  that 


It  follows  directly  from  the  recursion  formula  that 


(62) 


(63) 


(64) 


♦j  "  qJ  *  (65) 

The  three  Poisson  equations  for  the  velocity  components  are  solved  by 
sweeping  through  the  computational  domain  applying  the  recursion  formula  given 
by  Equation  (58).  The  manner  in  which  the  solution  of  the  Poisson-type  equa¬ 
tions  is  coupled  to  the  solution  of  the  transport-type  equations  is  discussed 
in  Section  3.4. 


3.2  Discretization  of  the  Trans port -Type  Equations 

As  given  in  Section  2,  the  transport  equations  for  the  vorticity  com¬ 
ponents,  the  turbulent  kinetic  energy,  and  the  dissipation  of  turbulent 
kinetic  energy  can  be  written  in  the  form 
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Vi  72  +  Vi  72  +  Vi  72  +  (V2  +  Vi}  ~= 

Y  3x  Y  dy  y  3z  r  y  3x 

+  <“*S2  +  W  ^  +  (V2  +  S*V  ^  "  %  ’ 


where  d>  *  Q  ,  fl  ,  G  ,  k,  or  e. 

T  x  y  z 

Equation  (66)  cannot  be  solved  for  an  arbitrary  Reynolds  number  using  the 
conventional  central-difference  approximations  to  the  derivatives.  The  reason 
for  this  difficulty  is  that  the  coefficients  of  the  convective  terms  3$/3x, 
34»/3y *  and  3<|>/3z  contain  the  Reynolds  number  as  a  multiplicative  factor.  Con¬ 
sequently,  with  the  standard  central-difference  algorithm,  the  discretized 
system  of  equations  is  diagonally  dominant  for  only  a  limited  range  in 
Reynolds  number.  Diagonal  dominance  is  necessary  to  obtain  convergence  in  the 
iterative  solution  of  the  discretized  system  of  equations. 

In  the  present  work  the  transport-type  equations  are  discretized  using 
Agarwal's  third-order-accurate  upwind-difference  scheme,  References  11  and 
12.  In  the  application  of  this  method.  Equation  (66)  is  rewritten  in  the  form 


(%a2  +  tyV 


V2 +  W 


(<*.Yi  +  fi.Y  1 ) 
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! 

i 

i 
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1 

'< 
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r  - 

fc‘ 
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The  basis  of  the  upwind-difference  algorithm  can  be  illustrated  by  considering 
the  term  3<j>/3x  in  Equation  (67).  At  an  interior  point  (i,j,k)  in  the  nodal 
network,  this  terra  is  evaluated  in  terras  of  the  adjacent  points  in  the  x 
direction  with  the  following  truncated  Taylor-series  representation  and  stand¬ 
ard  central-difference  approximation  to  the  first  derivative: 


dip 

3x 


i.J.k 


^1+1, j,k  ~  ^l-l.j.k  3^ 

2h  "6  3x3 


i.j.k 


iL 

"  51  3x5 


(74) 


i.j.k 


3  -3 

In  the  third-order-accurate  upwind-difference  scheme,  the  derivative  3  tp/dx 
is  retained  and  expressed  in  terms  of  one-sided  difference  approximations 
depending  on  the  sign  of  T : 


3x3 


324> 

3x2 


_  d^± 

~  3x2 


i-1 


if  rx  >  o. 


(75) 


334> 

3x3 


3 2_± 

3x2 


_ 

3x2 


i+1 


if  rx  <  o. 


(76) 


Using  Equations  (75)  and  (76),  it  follows  that 


-  r 


3$ 
1  3x 


„  /2*i+l  +  +  ♦i-2\  ,,  „ 

"  "  rl l  6h  f  '  1 1 


>  0 


(77) 


and 


-  r 


34> 
1  3x 


i+2  +  6*i+i  "  3*i  ”  2*i-l\  _ 

"  -  ri  [ - 6h - )  *  if  FI  <  0  * 


(78) 
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Analogous  formulas  apply  to  the  terms  -  (3$/3y)  and  -  (3<(>/3z)  appearing 

2  -2  2  ^  -2  2  -2 J 

in  Equation  (67).  The  terms  3  <|)/3x  ,  3  $/3y  ,  and  3  $/3z  and  the  derivatives 
appearing  in  the  coefficients  and  source  term  are  discretized  using  the  stand¬ 
ard  central-difference  approximations,  Equations  (38)  through  (43). 

With  this  approach,  the  following  is  the  finite-difference  form  of  the 
transport-type  equation: 

Tl*i+2,j,k  +  Vi+l,j,k  +  T3*i,j,k  +  Vi-l,j,k  +  T5*i-2,j,k 

+  Vi,j+2,k  +  Vi,j+l,k  +  Vi,j-l,k  +  Vi,j-2,k 

+  T10^i, j,k+2  +  Tll*i,j,k+1  +  T12*i,j,k-1  +  T13*i,j,k-2 


The  coefficients  Tj  through  Tjj  for  the  two  cases  >  0,  r2  >  0,  T-j  >  0  and 
^  r2  ^  r3  ^  ®  are  defined  in  Table  1  along  with  the  generalized  form 
of  the  coefficients  used  in  the  calculations.  When  Equation  (79)  is  applied 
at  the  second  node  from  a  boundary  plane,  the  value  of  $  at  the  boundary  is 
determined  from  the  boundary  condition  and  is  transferred  to  the  right  side  of 
the  equation.  When  Equation  (79)  is  applied  at  the  interior  node  adjacent  to 
the  boundary  plane,  the  value  of  $  at  the  boundary  is  determined  again  from 
the  boundary  condition,  and  the  value  of  $  at  the  external  node  adjacent  to 
the  boundary  is  determined  either  from  an  image  condition  or  extrapolation 
using  the  techniques  discussed  in  Section  3.3;  the  known  values  are  trans¬ 
ferred  to  the  right  side  of  the  equation.  To  solve  for  <j>  along  a  line  of 
variable  j,  Equation  (79)  is  written  in  the  matrix  form  as  follows: 
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TABLE  1.  COEFFICIENTS  IN  THE  DISCRETIZED  FORM  OF  THE  TRANSPORT  EQUATION. 


r,>o,  r2>o,  r3>o 

r,<o,  r2<o,  r3<o 

Generalized  Form 

0 

-  5  1  r>  1 

n  <ri  -  i  r>  i) 

A,  -  5  1  r,  1 

A,  +  h  |  r,  | 

a,  +  j  (1  r,  |  -  2r,) 

-(2A,  +  2A2  +  2Aj 

+  5  1  ri  1  +  \  1  r2 1 

—  (2  A)  +  2Aj  +  2Aj 

+  j  |r,|  +  \  1  r2| 

—  (2Aj  +  2A2  +  2Aj 

+  3  1 r, !  +  5  1  r2 1 

+  5  1  r3 1) 

+ 1  in  i) 

+  2  1  r3 1) 

a,  +  h  |  r,  | 

a,  -  j  1  r,  | 

A,  +  J  (|  r,  1  +  2T,) 

-  s  1  r«  1 

0 

-  n  (r,  +  1  r,  1) 

0 

-ghir2, 

n  (r2  -  1  r2 1) 

A2  -  5  1  r2 1 

a2  +  h  |  r2  | 

a2  +  5  (1  r2 1  -  2r2> 

a2  +  h  |  r2 1 

a2  -  5  1  r2 1 

a2  +  J  (i  r2 1  +  21^ 

-5'p*' 

0 

-  j5  <r2  +  1  r2 1) 

0 

-  6  i  r3  1 

n  <r3  -  1  r3 1) 

Aj  -  5  1  r3 1 

a3  +  h  1  r3  1 

a3  +  5  (1  r3 1  -  2r3> 

Aj  +  h  1  r3  1 

Aj  -  J  1  r3  I 

A,  +  J  (|  r3 1  +  2r3> 

-  S  1  r3  1 

0 

-  n  (p3  + 1  r3 1) 
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(80) 

Equation  (80)  represents  a  system  of  simultaneous  linear  equations  with  a 
coefficient  matrix  of  band  structure  having  a  total  band  width  of  five.  The 
system  is  solved  using  SUBROUTINE  GELB  from  the  McDonnell  Douglas  Automation 
Company  Scientific  Subroutine  Package,  Reference  13,  a  routine  which  performs 
Gaussian  elimination  with  column  pivoting  to  solve  the  matrix  problem. 

The  solution  sequence  of  the  transport-type  equations  is  discussed  in 
Section  3.4. 

3.3  Discretization  of  the  Boundary  Conditions 


As  discussed  in  Section  2.3,  boundary  conditions  imposed  in  the  flowfield 
analysis  are  of  either  the  Dirichlet  or  Neumann  type.  The  discretized  forms 
of  these  boundary  conditions  are  discussed  in  this  section. 
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On  an  inflow  plane  where  the  velocity  and  turbulence  property  profiles 
are  specified  analytically,  the  boundary  conditions  given  by  Equation  (28)  are 
expressed  as  known  functions  of  the  discretized  coordinates, 

ui  J  k  =  u^xi»  zk)>  vi  j  k  "  v(xi»  wi  j  k  ’  w(xi»  zic^* 

ki  j  k  -  k^,  h,  z^),  and  ei  j  k  ■  eCx^  K,  zfe).  (81) 

The  vorticity  profiles  are  computed  by  differentiating  the  analytic  velocity 
profiles  according  to  Equations  (29)  through  (31),  and  the  resulting  expres¬ 
sions  are  given  as  functions  of  the  discretized  coordinates, 


ft  -  ft  (x  ,  h,  z  ),  ft  -  ft  (x  ,  h,  z  )  , 

xi,J,k  x  1  K  yi,J,k  y  1  K 

and  ft  -  ft  (x , ,  ft,  z.  ).  (£ 

zi.j.k  ’  1  k 

On  solid  surfaces  the  velocity  components  and  turbulent  kinetic  energy 
are  expressed  in  finite-difference  form;  for  example,  on  the  ground  plane 


Ui,2,k  “  Vi,2,k  "  Wi,2,k  "  ki,2,k  "  °* 


On  Che  ground  plane  the  turbulent  dissipation  follows  directly  from  Equation 
(33)  using  a  one-sided  difference  formula  to  represent  the  second  derivative 
of  the  turbulent  kinetic  energy. 


The  vorticity  components  on  a  solid  surface  are  discretized  using  the 
approach  developed  by  Woods,  Reference  14.  Consider  ft^  evaluated  on  the 
surface  y  ■  0.  Since  v  =  0  on  the  ground  plane,  Equation  (29)  reduces  to 


2  2 

It  follows  that  Sft^/dy  ■  3  w/3y  ,  or  in  terms  of  the  transformed  y  coordinate, 


9  5L.B2^+6  i«. 

3y  3y  3y 


(86) 


Using  a  Taylor  series,  w  at  the  node  adjacent  to  the  wall  is  represented  by 
the  following  expansion: 


W  m  V  + 

i,3,k  i,2,k 


3w 

3y 


i,2,k 


2»? 


(87) 


i,2,k 


When  Equations  (83),  (84),  and  (85)  are  introduced,  Equation  (87)  becomes 


w 


i,2,k 


.2  3ft 


i,3,k  g 


1,2 


2  6 


1,2  3y 


S2,2  1,2  fl« 


i,2,k 


(88) 


i,2,k 


'1,2 


A  one-sided  difference  formula  is  used  to  represent  the  vorticity  derivative. 


3ft  flx.  .  .  flx.  .  . 

_ x  i,3,k _ i,2,k 

3y  h 

and  Equation  (88)  is  solved  for  the  wall  vorticity: 


(89) 


2  ei,2  Wi,3,k  ‘  01,2  nx 


l,3,k 


i,2,k 


01,2  h  ”  02,2  b2 


(90) 


On  the  ground  plane 


ft  -  0 

yi,2,k 


(91) 
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and 


2  3 l ,2  Ui,3,k  +  0l,2  nz 


i,3,k 


"i»2,k 


^2,2  k  -  ^1,2  h 


On  outflow  and  symmetry  planes  where  the  normal  gradients  vanish,  the 
derivatives  expressed  by  Equations  (34)  through  (36)  are  represented  in 
discrete  form  using  a  one-sided  finite-difference  formula.  For  example, 
consider  the  condition  3$/3x  =  0  evaluated  on  the  plane  x  =  0  (i  =  2): 


*3,j,k  ~  H.j.k  ~  3<^2,j,k  _ 

2h  _  U» 


2,j,k 


which  results  in  the  boundary  value  of  $  being  given  by 


*2,j,k  “  3  ^3, j,k  "  3  *4,j,k  * 


In  the  solution  of  the  transport-type  equations,  values  of  the  flow 
properties  are  required  on  planes  external  to  the  boundaries  of  the  compu¬ 
tational  domain  and  spaced  a  distance  h  from  the  latter.  If  the  gradient  of 
the  flow  property  vanishes  on  the  computational  boundary,  the  external  plane 
is  an  image  plane  and  the  value  of  the  property  is  determined  by  the  vanishing 
gradient.  For  example,  at  the  surface  x  ■  0  (i  *  2), 


*l,j,k  ~  »3,j,k 


2*  j,k 


*3, j,k  * 

If  the  gradient  does  not  vanish,  the  flow  property  on  the  external  plane  is 
evaluated  from  a  quadratic  extrapolation  formula.  For  example, 


*l,j,k  "  3*2,j,k  "  3*3,j,k  +  *4,j,k  * 
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In  discretization  of  the  boundary  conditions,  special  attention  had  to  be 
given  to  evaluating  velocity  components  normal  to  the  outflow  boundaries. 
Through  the  use  of  the  coordinate  transformations,  the  boundary  conditions  on 
the  outflow  surfaces  are  specified  in  the  far  field.  The  level  of  the  normal 
velocity  component  is  not  specified;  rather,  its  normal  derivative  is  set  to 
zero,  and  the  level  of  the  velocity  is  computed  as  part  of  the  solution.  It 
was  found  that  following  the  achievement  of  a  convergent  flowfield,  an  inte¬ 
gral  mass  balance  performed  on  the  computational  domain  failed  to  satisfy  con¬ 
tinuity  within  an  acceptable  tolerance.  Consequently,  a  technique  was  devised 
to  modify  the  normal  velocity  component  on  each  outflow  plane  in  the  manner 
described  below. 


The  mass  flow  rate  across  each  boundary  of  the  solution  domain  is  com¬ 
puted  using  trapezoidal  integration.  For  example,  on  the  plane  x  «  L,  the 
normalized  mass  flow  rate  is  given  by 


m_ 

x=L 


lu|  cos  8 
BiTi 


dzdy 


(u  >  0,  cos  9  ■  1;  u  <  0,  cos  9  *  -  1).  (98) 

A  ratio  is  then  formed  of  the  net  flow  into  the  computational  domain  to  the 
net  flow  from  the  computational  domain,  and  the  normal  velocity  component 
computed  from  the  formula  given  by  Equation  (94)  is  multiplied  by  the  ratio 
every  10  iterations  in  the  solution  cycle.  As  the  flow-variable  residuals 
decline,  the  mass  flow  rate  approaches  unity  and  the  net  mass  flow  rate 
reduces  to  10-^. 

Specific  boundary  conditions  for  the  flow  configurations  considered  are 
defined  in  Section  4. 


3.4  Solution  of  the  Coupled  System  of  Equations 

The  sequence  in  which  the  solution  of  the  Poisson-  and  transport-type 
equations  is  carried  out  is  illustrated  in  Figure  4.  Initial  distributions 
are  established  for  all  flow  properties.  The  Poisson  equations  are  solved  for 
the  velocity  components  using  the  tridiagonal  algorithm  and  alternating- 
direction  line  iteration  in  which  the  values  of  the  velocity  components  are 
updated  according  to  the  formula 


29 


I 


4> 


n+1 

i.  j»k 


♦J'1.  .  +  R 

yi,J,k 


-  <(> 


n-1 

i,  j» 


(99) 


where  n  is  the  iteration  counter  and  R  is  the  relaxation  factor,  which  has  a 
value  of  1.2  for  the  Poisson  equations.  The  transport-type  equations  are 
solved  for  the  vorticity  components  and  turbulence  quantities  using  the  penta- 
diagonal  algorithm  (GELB)  and  alternating-direction  line  iteration  in  which 
the  properties  are  updated  according  to  Equation  (99)  with  R  ranging  from  0.2 
to  0.8  depending  upon  the  jet  conf iguration.  As  soon  as  a  flow  variable  has 
been  recomputed,  the  associated  boundary  conditions  are  re-evaluated.  Con¬ 
vergence  is  considered  to  have  been  achieved  when  the  root-mean-square  resi¬ 
dual  for  each  variable. 


r 


(100) 


has  reached  a  specified  level  (10-^  to  10“^  depending  upon  the  magnitude  of 
the  unknown). 


The  flowfield  calculations  in  the  present  work  were  done  using  the  CDC 
7600  computer  of  the  Systems  Technology  Program  -  System  Simulation  Center 
operated  by  the  McDonnell  Douglas  Astronautics  Company,  Huntington  Beach,  for 
the  Department  of  the  Army. 

Flowfields  were  computed  for  isolated  and  interacting  jets.  The  maximum 
Reynolds  number  for  which  convergent  solutions  could  be  achieved  was  on  the 
order  of  10^  using  the  finest  finite-difference  grid  permitted  by  the  avail¬ 
able  storage  capacity  of  the  CDC  7600.  Each  of  these  solutions  required 
approximately  1  h  of  machine  computation  time.  Definitions  of  the  flow  con¬ 
figurations  and  comparisons  of  the  calculated  properties  with  data  are 
presented  in  the  next  section. 
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4.  THE  COMPUTED  FLOWFIELDS 


Under  the  present  contract,  turbulent  flowfields  were  calculated  for 
Isolated  and  interacting  jets  in  ground  effect.  For  the  single  jet,  both 
normal  and  inclined  impingement  were  considered,  and  for  the  two  jets,  both 
equal  and  unequal  nozzle  diameters  were  considered.  In  this  section,  calcu¬ 
lated  flow  variables  are  presented  for  these  configurations  (as  well  as  for  a 
laminar-flow  test  problem),  and  comparisons  are  made  between  computed  and 
measured  properties. 

4.1  Laminar  Channel  Flow 

In  order  to  establish  the  accuracy  of  the  numerical  algorithm  for  solu¬ 
tion  of  the  Navier-Stokes  equations,  laminar  flow  in  a  channel  was  computed  as 
one  of  the  standard  problems  considered  in  evaluation  of  any  new  finite- 
difference  method.  A  quarter-section  of  the  channel  is  shown  in  Figure  5(a), 
where  the  half-width  is  denoted  by  W,  the  half-height  by  H,  and  the  length  by 
L.  Since  the  entering  flow  is  taken  to  be  symmetric  with  respect  to  the 
channel  midplanes  y  •  0  and  z  *  0,  only  a  quarter  section  of  the  passage  need 
be  computed.  In  this  problem,  the  aspect  ratio  of  the  channel  is  unity  (W  * 

H)  and  L  is  equal  to  20W.  Normalizing  by  the  channel  width  2W,  the  dimensions 
of  the  physical  solution  domain  in  the  x,  y,  and  z  directions  are  10.0,  0.5, 
and  0.5,  respectively. 

For  this  problem,  a  coordinate  transformation  is  used  in  the  x  direction 
only,  where  the  computational  coordinate  x  (normalized  by  the  channel  width) 
is  related  to  the  physical  x  by  the  following  expression: 


Figure  5.  Velocity  profiles  for  laminar  channel  flow  (Re  =  75). 
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Values  of  x  for  uniform  increments  of  x  are  tabulated  below. 


X 

X 

X 

X 

0. 

0. 

1.1 

1.0000 

0.1 

0.0476 

1.2 

1 . 2000 

0.2 

0.1000 

1.3 

1.4444 

0.3 

0.1579 

1.4 

1.7500 

0.4 

0.2222 

1.5 

2.1429 

0.5 

0.2941 

1.6 

2.6667 

0.6 

0.3750 

1.7 

3.4000 

0.7 

0.4667 

1.8 

4.5000 

0.8 

0.5714 

1.9 

6.3333 

0.9 

0.6923 

2.0 

10.0000 

1.0 

0.8333 

Thus,  the  dimensions  of  the  computational  domain  in  the  x,  y,  and  z 
directions,  respectively,  are  2.0  x  0.5  x  0.5. 

The  boundary  conditions  for  the  velocity  components  are  given  in  Figure 
5(a).  At  the  entrance  plane  to  the  channel,  it  is  assumed  that  the  flow  is  in 
the  x  direction  only.  Symmetry  conditions  are  imposed  on  the  planes  y  ■  0  and 
z  =»  0,  and  the  no-slip,  impermeable-wall  constraints  are  enforced  at  the  solid 
surfaces  x  =  W  and  y  =  H.  At  the  outflow  plane,  the  velocity  components  are 
taken  to  be  invariant  with  x  (fully  developed  flow). 

The  three-dimensional  laminar  channel  flow  was  computed  with  Re  ■  75 
(based  on  the  channel  width  and  entrance-plane  velocity)  using  42  grid  points 
in  the  x  direction,  12  in  the  y  direction,  and  12  in  the  z  direction.  The 
centerline  velocity  variation  with  the  axial  coordinate  is  shown  in  Figure 
5(b),  and  the  streamwise  velocity  profiles  on  the  midplane  z  *  0  for  selected 
x  stations  are  shown  in  Figure  5(c).  The  slight  difference  between  the 
numerical  value  of  the  fully  developed  centerline  velocity  and  that  determined 
analytically  is  attributed  to  the  relative  coarseness  of  the  mesh,  as  is  the 
small  disparity  between  the  numerical  streamwise  velocity  profile  for  x  =  0.57 
and  the  data  of  Goldstein  and  Kreid,  Reference  15.  Two  contour  plots  of  the 
x-component  of  velocity  are  shown  in  Figure  5(d)  which  demonstrate  the 
acceleration  of  the  central  core  as  the  boundary  layers  develop  on  the 
confining  walls. 

4.2  Turbulent  Isolated  Jet-Impingement  Flow 


The  first  turbulent  jet-impingement  flow  considered  is  that  for  an  iso¬ 
lated  nozzle  with  a  free  upper  boundary  and  its  axis  perpendicular  to  the 
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ground,  as  illustrated  in  Figure  1(a)  in  planar  view.  The  flowfield  is  estab¬ 
lished  with  specification  of  the  nozzle  height  above  ground,  the  velocity  and 
turbulence  profiles  at  the  nozzle  exit  plane,  and  the  Reynolds  number  based  on 
the  nozzle  exit  conditions. 

In  order  to  take  advantage  of  symmetry,  only  a  quarter  of  the  jet  flow- 
field  is  computed,  as  illustrated  in  Figure  6(a).  The  top  of  the  computa¬ 
tional  domain  is  not  at  the  nozzle  exit  plane  but  rather  is  located  a 
distance  K  above  the  ground.  The  solution  domain  extends  a  distance  L  in  the 
x  direction  and  a  distance  W  in  the  z  direction.  In  this  problem  the  width 
and  depth  are  taken  to  be  equal  (L  =*  W),  and  the  dimensions  of  the  physical 
domain  normalized  by  the  jet  nozzle  diameter  are  chosen  to  be  L  *  2.0,  K  » 

1.5,  and  W  *  2.0. 

Coordinate  transformations  are  used  in  both  the  x  and  z  directions,  where 
the  computational  coordinates  x  and  z  (normalized  by  the  nozzle  diameter)  are 
expressed  in  terras  of  the  physical  coordinates  through  the  following 
relations: 


x  =,  2-25x 

1  +  x 


(102) 


-  2.25z 

Z  3  ■  • 

1  +  Z 


Values  of  x  and  z  for  uniform  increments 


x,  z 

0. 

0.1 

0.2 

0.3 

0.4 

0.5 

0.6 

0.7 

0.8 

0.9 

1.0 

1.1 

1.2 

1.3 

1.4 

1.5 


(103) 


of  x  and  z  are  listed  below. 


x,  z 

0. 

0.0465 

0.0976 

0.1538 

0.2162 

0.2857 

0.3636 

0.4516 

0.5517 

0.6667 

0.8000 

0.9565 

1.1429 

1.3684 

1.6471 

2.0000 
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The  dimensions  of  Che  computational  domain  are  1.5  x  1.5  x  1.5  in  the 
x,  y  and  z  directions,  respectively. 


The  velocity  profile  at  the  inflow  plane  is  assumed  to  be  fully  developed 
and  given  by  the  relation 


v(x,  h,  z)  =  - 


1 


1  +  Kr 


2  ’ 


(104) 


where  K  is  a  constant  for  a  fixed  value  of  ft  and  r  is  the  distance  along  the 
perpendicular  from  points  within  the  inflow  plane  to  the  jet  centerline: 


K  = 


0.414 


(105) 


(5/D)' 


2  2 
x  +  z  . 


(106) 


In  Equation  (105)  6  is  the  jet  half-radius  at  a  distance  (H  -  h)  below  the 
nozzle  exit  plane  and  is  computed  from  the  empirical  relations  of  Giralt, 
Chia,  and  Trass,  Reference  16.  Also  on  the  inflow  plane,  v  *  w  =  0,  and  k 
and  e  follow  from  the  assumed  profiles  given  below. 


k(x,  h,  z)  =  -  0.04v(x,  h,  z) 


e(x,  h,  z)  ■  5k(x,  h,  z) 


3/2 


(107) 

(108) 


On  the  impingement  surface  the  velocity  components  and  the  turbulent  kinetic 
energy  are  zero,  and  the  dissipation  is  determined  from  the  reduced  form  of 
the  turbulent-kinetic-energy  transport  equation  at  the  wall.  The  boundary 
conditions  follow  from  symmetry  conditions  for  x  =  0  and  z  =  0  and  from  the 
assumption  of  developed  flow  for  x  *  L  and  z  =  W. 

The  impinging- jet  flowfield  was  computed  for  Re  =  100  and  H  *  7.5  (h  * 
1.5)  using  17  grid  points  in  each  of  the  coordinate  directions.  Figure  6(b) 
illustrates  the  ground-plane  pressure  variation  normalized  by  the  stagnation- 
point  pressure  p„,  where  p  is  the  minimum  pressure  on  the  surface.  The 
computed  profile,  which  has  been  reflected  about  the  symmetry  plane  x  =  0,  is 
compared  with  the  experimental  data  of  Snedeker  and  Donaldson,  Reference  17, 
for  the  same  value  of  H  but  for  a  Reynolds  number  on  the  order  of  100  000. 


Figure  6.  Flow  properties  for  an  isolated  jet  with  normal  impingement. 
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For  the  same  configuration,  Figure  7(a),  the  ground-plane  pressure  and 
jet  centerline  velocity  as  a  function  of  distance  from  the  stagnation  point 
normalized  by  the  jet  half-radius  at  the  beginning  of  the  impingement  region 
are  illustrated  in  Figures  7(b)  and  7(c),  respectively,  and  are  compared  with 
the  experimental  data  of  Giralt,  Chia,  and  Trass,  Reference  15.  The  length 
scale  6  is  clearly  useful  since  it  collapses  the  pressure  and  velocity 
profiles  to  a  single  curve  for  a  range  of  H.  The  ground-plane  pressure 
profile  for  z  =  0  computed  in  the  present  work  with  Re  =  100  compares 
reasonably  well  with  the  experimental  results,  which  are  not  strongly 
dependent  on  Reynolds  number  since  the  surface  pressure  in  the  impingement 
region  is  governed  by  the  dynamics  of  the  flow.  For  this  reason,  the  inviscid 
solution  of  Scholtz  and  Trass,  Reference  18,  which  is  also  shown  in  Figure 
7(b),  agrees  well  with  the  data.  The  dimensionless  centerline  velocity  decay 
of  the  impinging  jet  is  plotted  in  Figure  7(c)  as  a  function  of  the  distance 
from  the  stagnation  point  normalized  by  the  half-radius  of  the  jet  at  the 
beginning  of  the  impingement  region.  The  present  solution  for  Re  =  100 
follows  the  experimental  centerline  velocity  variation  obtained  for  Reynolds 
numbers  of  30  000  to  80  000.  However,  the  velocity  decay  computed  by  Scholtz 
and  Trass,  Reference  18,  using  inviscid  theory  overpredicts  the  magnitude  of 
v.  The  reason  for  this  discrepancy  is  that  for  H/D  >  7.78,  as  is  the  case  in 
the  measurements,  the  impinging- jet  velocity  decay  and  spreading  are 
appreciably  affected  by  turbulent  free-jet  entrainment  and  viscous  effects. 

For  values  of  H/D  <  7.78,  the  accuracy  of  the  inviscid  theory  should  improve, 
however,  since  the  free-jet  entrainment  is  small  in  such  cases  and  the  jet 
decay  and  spreading  during  the  deflection  are  primarily  governed  by  the 
presence  of  the  wall. 

A  more  complex  Isolated  jet-impingement  flow  is  formed  when  the  nozzle 
centerline  is  inclined  at  an  angle  9  with  respect  to  the  ground  plane.  The 
computational  domain  for  this  configuration  is  illustrated  in  Figure  8(a).  In 
this  case  there  is  only  a  single  symmetry  plane,  and  half  the  jet  flowfield 
must  be  computed.  The  inflow  plane  is  located  at  y  ■  h;  the  outflow  planes 
are  at  x  *  0,  x  =  L,  and  z  =  W;  and  the  symmetry  plane  is  at  z  =  0.  For  this 
problem  L  =  4,  h  =  1.5,  and  W  =  2. 
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The  following  coordinate  transformations  are  used  in  the  x  and  z 
directions: 


0.75x 

X  *  -r - 

3  -  X 


0  <  x  <  1.5 


(109) 


"  „  3.75x  -  6 
x  -  1 


1.5  <  x  <  3.0 


(110) 


2.25z 
1  +  z  * 


(111) 


Tabulated  values  of  x  and  z  for  uniform  increments  in  x  and  z  are  given  below. 


X 

X 

X 

X 

0. 

0. 

1.6 

2.0465 

0.1 

0.3529 

1.7 

2.0976 

0.2 

0.6316 

1.8 

2.1538 

0.3 

0.8571 

1.9 

2.2162 

0.4 

1.0435 

2.0 

2.2857 

0.5 

1.2000 

2.1 

2.3636 

0.6 

1.3333 

2.2 

2.4516 

0.7 

1.4483 

2.3 

2.5517 

0.8 

1.5484 

2.4 

2.6667 

0.9 

1.6364 

2.5 

2.8000 

1.0 

1.7143 

2.6 

2.9565 

1.1 

1.7838 

2.7 

3.1429 

1.2 

1.8462 

2.8 

3.3684 

1.3 

1.9024 

2.9 

3.6471 

1.4 

1.9535 

3.0 

4.0000 

1.5 

2.0000 

_z 

z 

0. 

0. 

0.1 

0.0465 

0.2 

0.0976 

0.3 

0.1538 

0.4 

0.2162 

0.5 

0.2857 

0.6 

0.3636 

0.7 

0.4516 

0.8 

0.5517 

0.9 

0.6667 

1.0 

0.8000 

1.1 

0.9565 

1.2 

1.1429 

1.3 

1.3684 

1.4 

1.6471 

1.5 

2.0000 

*1 
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The  dimensions  of  the  computational  domain  in  the  x,  y,  and  z  directions, 
respectively,  are  3.0  x  1.5  x  1.5. 

The  velocity  profile  at  the  inflow  plane  is  again  assumed  to  be  fully 
developed  and  is  given  by  the  relation 

V(r)  - - l—^r  .  (112) 

1  +  Kr 

In  Equation  (112),  K  is  the  constant  defined  by  Equation  (105)  and  r  is  the 
distance  along  the  perpendicular  from  points  within  the  inflow  plane  to  the 
jet  centerline,  which  in  this  case  is  given  by 

r^  *  [(x  -  L/2)  sin  6  -  y  cos  0]^  +  z^  .  (113) 

The  boundary  values  for  the  velocity  components  are  given  by  the  relations: 


u(x,  h. 

z)  *  V(r)  cos  0 

(114) 

v(x,  h. 

z)  ■  V(r)  sin  0 

(115) 

w(x,  h, 

z)  =  0. 

(116) 

The  turbulent  kinetic  energy  and  dissipation  are  computed  from  Equations  (107) 
and  (108).  The  remaining  boundary  conditions  on  the  impingement  surface  and 
the  three  outflow  planes  are  defined  in  Figure  8(a). 

The  flowfield  for  this  configuration  was  computed  for  Re  ■  100  and  H  ■ 

7.5  (h  =  1.5)  using  32  grid  points  in  the  x  direction,  17  in  the  y  direction, 
and  17  in  the  z  direction.  The  computed  impingement-plane  pressure  along  the 
line  z  *  0  is  shown  in  Figure  8(b)  for  9  =  90°,  75°,  and  60°  and  compared  with 
the  data  of  Snedeker  and  Donaldson,  Reference  17,  for  the  case  of  inclined 
impingement.  In  this  plot,  p  is  the  minimum  surface  pressure;  p_  _  is  the 
stagnation-point  pressure  for  normal  impingement;  and  x'  *  x  -  L/2.  The  com¬ 
puted  variations  of  the  stagnation-point  pressure  and  displacement  as  func¬ 
tions  of  the  jet  vector  angle  agree  reasonably  well  with  the  data  as  shown  in 
Figures  8(c)  and  8(d).  For  the  same  inclined-jet  configuration  shown  in 
Figure  9(a),  Figures  9(b)  through  9(e)  illustrate  the  y  component  of  velocity 
and  the  x  component  of  velocity  along  the  plane  z  •  0  for  9  ■  75°  and  9  ■  60°. 
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(b)  y -component  of  velocity  for  i  »  0,  6  «  75°  (c)  y -component  of  velocity  for  z  =  0,  6  =  60° 


t.pjt-een-i.' 


Figure  9.  Flow  properties  for  ■■  Isolated  jet  with  IncHned  impingement. 


A  more  complicated  jet-impingement  flow  exists  when  the  jet  discharges 
through  a  blocking  plate  whic'i  is  parallel  and  in  proximity  to  the  ground.  In 
this  case  the  flow  separates  on  the  upper  surface,  and  a  recirculating  region 
must  be  treated  in  the  analysis.  The  computational  domain  is  illustrated  in 
figure  10(a),  where  a  quarter  of  the  jet  is  computed  for  normal  impingement. 
The  spacing  between  the  plates  is  denoted  by  H,  and  the  plate  dimensions  are 
WxL,  where  all  lengths  are  normalized  by  the  diameter  of  the  entering  jet. 
Taking  L  =  W,  the  dimensions  of  the  physical  solution  domain  in  the  x,  y,  and 
z  directions  are  5.0  x  1.5  x  5.0,  respectively. 

For  this  problem,  coordinate  transformations  are  used  in  both  the  x  and  z 
directions,  where  the  computational  coordinates  x  and  z  are  expressed  in 
terras  of  the  physical  coordinates  x  and  z  by  the  following  relations: 


x 


1 .8x 
1  +  x 


(117) 


-  _  1.8z 

1  +  z  * 


(118) 


Tabular  values  of  x  and  z  for  uniform  increments  of  the  transformed 
coordinates  are  listed  below. 


X,  z 

x,  z 

0. 

0. 

0.1 

0.0588 

0.2 

0.1250 

0.3 

0.2000 

0.4 

0.2857 

0.5 

0.3846 

0.6 

0.5000 

0.7 

0.6364 

0.8 

0.8000 

0.9 

1.0000 

1.0 

1.2500 

1.1 

1.5714 

1.2 

2.0000 

1.3 

2.6000 

1.4 

3.5000 

1.5 

5.0000 

The  dimensions  of  the  computational  domain  are  1.5  x  1.5  x  1.5 


(•)  Computational  domain 


H:u  =  0,  v  -  v  (x.H.  z).w  =  0,  k  =  k  (x.  H.  z). 
e  =  €  (x,  H,  z) 


<,pj  t-aan-M 


Figure  10.  Flow  properties  for  in  isolated  Jet  with  normal  Impingement  in  the  presence  of  a  blocking  plate. 
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The  boundary  conditions  for  the  velocity  and  turbulence  quantities  are 
defined  in  Figure  10(a).  The  entering  jet  velocity  profile  at  the  opening  in 
the  top  plate  is  specified  according  to  the  following  relation: 


v(x,  H,  z)  -  -  1.0 


r  <  r 


(119) 


v(x,  H,  z) 


(r  -  re)2  -  (r£  -  rfa)2 
(re  -  V2 


r  <  r  <  r.  , 
e  b 


(120) 


where  r  =  (x^  +  z2)l/2^  r  _  0.5^  an(j  was  chosen  to  be  0.4  in  the  present 
calculations.  For  r  >  0.5,  v(x,  H,  z)  =  0;  u  and  w  are  zero  along  the  entire 
upper  surface  of  the  computational  domain.  Within  the  entering  jet,  k  and  e 
are  computed  using  Equations  (107)  and  (108);  outside  the  jet  along  the  upper 
plate,  the  turbulent  kinetic  energy  vanishes,  and  the  dissipation  is  computed 
from  the  wall  form  of  the  k  transport  equation, 


_1  a2k 
Re  2  * 

ay 


(121) 


Symmetry  conditions  are  applied  on  the  planes  x  »  0  and  z  *  0,  and  the  deriva¬ 
tives  of  all  flow  properties  with  respect  to  the  coordinate  normal  to  the  sur¬ 
face  are  set  equal  to  zero  on  the  outflow  planes  x  =  L  and  z  =  W.  On  the 
impingement  surface,  the  three  velocity  components  and  the  turbulent  kinetic 
energy  are  set  to  zero,  and  the  dissipation  is  evaluated  using  Equation  (121). 

The  flowfield  was  calculated  for  the  three-dimensional  impinging  jet  in 
the  presence  of  a  blocking  plate  with  Re  »  100  and  using  17  grid  points  in 
each  coordinate  direction.  A  streamline  pattern  for  this  configuration  given 
in  the  computational  (not  physical)  solution  domain  is  shown  in  Figure  10(b), 
where  pathlines  are  traced  from  the  boundary  of  the  entering  jet.  The  pro¬ 
files  of  the  y  component  of  velocity  on  the  symmetry  plane  x  *  0  are  given  in 
Figure  10(c)  and  illustrate  the  decay  of  the  jet  velocity  as  the  ground  is 
approached.  The  computed  ground-plane  pressure  on  the  line  x  =*  0  is  shown  in 
Figure  10(d)  along  with  the  data  of  Bradbury,  Reference  19,  for  the  normal 
impingement  of  an  axisymmetrlc  jet  in  the  absence  of  a  blocking  surface. 
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4 . 3  Turbulent  Interacting  Jet -Impingement  Flow 

The  jet-impingement  flow  of  primary  interest  in  the  present  study  con¬ 
sists  of  two  jets  impinging  perpendicular  to  the  ground,  creating  opposing 
wall  jets  that  collide  to  form  a  fountain  upwash,  as  illustrated  in  Figure 
1(b).  Computations  were  made  for  this  configuration  with  jets  of  equal  and 
unequal  initial  diameters  discharging  through  a  blocking  plate.  The  flowfield 
is  established  with  specification  of  the  nozzle  height  above  ground,  the 
centerline  spacing,  the  diameter  ratio,  the  velocity  and  turbulence  profiles 
at  the  nozzle  exit  plane,  and  the  Reynolds  number  based  on  the  nozzle  exit 
conditions. 

When  the  jets  are  of  equal  diameter  and  equal  strength,  the  flowfield  is 
symmetrical  with  respect  to  the  midplane  of  the  fountain,  and  only  half  of  the 
flow  need  be  computed.  The  computational  domain  and  the  associated  boundary 
conditions  used  in  this  case  are  the  same  as  those  shown  in  Figure  10(a)  with 
L  =  S/2  and  the  boundary  condition  3u/3x  *  0  for  x  -  L  replaced  by  u  =  0.  The 
entering  jet  velocity  profile  is  given  by  Equations  (119)  and  (120). 

However,  the  coordinate  transformation  in  the  x  direction  used  for  this 
geometry  is  modified  from  that  used  for  the  isolated  jet  to  provide  the  proper 
grid  resolution  in  the  vicinity  of  the  fountain  symmetry  plane.  The  following 
transformation  is  used: 


x  = 


Ax 


Bx  +  C  ’ 


x  =  S  ,  +  S  „ 
cl  c2 


0  <  x  <  S  . 

cl 

D(Spl  +  Sp2  '  *> 

^pT  *  Sp2  ’ 


(122) 


S,<x<S,+S0, 
cl  cl  c2 


(123) 


where 


C  =  1 


(126) 


I 

! 


i 


D  = 


Sc2Spl(Scl  +  Sc2) 
■SclSp2(Spl  +  Sp2) 


E  =  Sc2Spl  ~  SclSP2 

'  Sc‘lSP2(Spl  +  S~PV 


F  =  1  . 


(127) 


(128) 


(129) 


In  Equations  (122)  through  (129),  x  =  Scj  denotes  the  point  in  the  com¬ 
putational  domain  where  a  coordinate  transformation  of  the  form  given  by 
Equation  (122)  originating  at  x  ■  0  is  matched  to  a  coordinate  transformation 
of  the  same  form  originating  at  the  fountain  midplane  (x  ■  Scj  4-  S^)  and 
extending  toward  the  junction  point.  The  constants  A  through  F  are  determined 
such  that  their  values  and  first  two  derivatives  of  the  two  transformations 
are  continuous  at  x  >  Scp  The  line  x  *  Scl  in  the  computational  plane 
corresponds  to  the  line  x  »  Sp^  in  the  physical  domain,  and  the  line  x  *  + 

SC2  in  the  computational  domain  corresponds  to  the  line  x  =  Sp^  +  Sp2  in  the 
physical  domain. 


For  a  twin 

-jet  configuration  with  a 

centerline 

spacing  of  6  nozzle 

diameters,  the 

following  values  were  used 

in  the 

transformation:  Scj  « 

Sc2  =  °*5»  Spl 
tabulated  below 

■  2.0,  and  Sp2  =  1.0.  The 

values 

of 

x  as  a  function  of 

X 

X 

X 

x 

0. 

0. 

0.1 

0.1017 

1.1 

1.3469 

0.2 

0.2069 

1.2 

1.5000 

0.3 

0.3158 

1.3 

1.6596 

0.4 

0.4286 

1.4 

1.8261 

0.5 

0.5454 

1.5 

2.0000 

0.6 

0.6667 

1.6 

2.1818 

0.7 

0.7924 

1.7 

2.3721 

0.8 

0.9231 

1.8 

2.5714 

0.9 

1.0588 

1.9 

2.7805 

1.0 

1.2000 

2.0 

3.0000 
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The  coordinate  transformation  given  by  Equation  (118)  is  used  in  the  z  direc¬ 
tion;  no  transformation  is  used  in  the  y  direction.  The  dimensions  of  the 
physical  domain  are  3.0  x  3.0  x  5.0  in  the  x,  y,  and  z  directions,  respec¬ 
tively,  and  the  dimensions  of  the  computational  domain  are  2.0  x  3.0  x  1.5. 

The  flowfield  was  computed  for  this  configuration  (S  =  6  and  H  =  3), 
which  is  presented  in  planar  view  in  Figure  11(a),  with  Re  =  200  and  22  grid 
points  in  the  x  direction,  32  in  the  y  direction,  and  17  in  the  z  direction. 
The  jet  centerline  variations  of  the  velocity  and  turbulent  kinetic  energy  are 
illustrated  in  Figure  11(b),  and  contour  plots  of  the  x  and  y  components  of 
velocity  are  shown  for  the  planes  z  ■  0  ana  z  =  0.5  in  Figures  11(c)  through 
11(f).  In  these  plots  x"*  ■  x  -  S/2.  The  ground-plane  pressure  distribution 
along  the  line  z  *  0  is  shown  in  Figure  12(a),  where  the  computed  profile  has 
been  reflected  about  the  fountain  symmetry  plane  and  compared  with  the  data  of 
Jenkins  and  Hill,  Reference  20,  for  the  same  configuration  in  the  absence  of  a 
blocking  plate.  The  calculations  provide  good  correlation  with  data  in  the 
impingement  zone  where  the  pressure  gradient  is  associated  with  the  deflection 
of  the  jet  by  the  ground  plane  and  is  not  significantly  influenced  by  viscous 
effects.  However,  in  the  fountain  region,  viscous  effects  determine  the  level 
of  the  pressure  recovery.  (In  an  inviscid  analysis  of  the  present  problem, 
the  pressure  would  be  unity  at  the  fountain  midline  x'  »  0,  y  »  0  since  it  is 
a  stagnation  line.)  For  the  relatively  low  Reynolds  number  in  the  calcula¬ 
tions,  diffusive  effects  are  large  and,  as  illustrated  in  Figure  12(b),  the 
opposing  wall  jets  rapidly  lose  momentum  such  that  the  computed  pressure 
increase  in  the  vicinity  of  the  fountain  is  smaller  than  the  measured  increase 
given  in  Reference  20,  which  corresponds  to  a  Reynolds  number  on  the  order  of 
100  00 0.  Calculations  were  first  made  with  Re  *  100,  and  there  was  no  appre¬ 
ciable  ground-plane  pressure  increase  in  the  fountain.  This  solution  was 
taken  to  initiate  the  solution  for  Re  -  200,  and  using  the  latter,  a  solution 
for  Re  =  300  was  attempted;  however,  a  convergent  solution  could  not  be 
achieved  for  the  highest  Reynolds  nunber. 

The  flowfield  for  equal-strength  jets  with  S/D  *  12  and  H/D  =  4  was  com¬ 
puted  using  the  x-coordinate  transformation  given  by  Equations  (122)  and  (123) 
with  Sc^  -  2.2,  Sc2  *  0.8,  Spl  ■  4.0,  and  Sp2  M  2.0.  The  values  of  the 
transformed  and  physical  coordinates  are  listed  below. 
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» 


X 

X 

0. 

0. 

0.1 

0.1468 

0.2 

0.2963 

0.3 

0.4486 

0.4 

0.6038 

0.5 

0.7619 

0.6 

0.9231 

0.7 

1.0874 

0.8 

1.2549 

0.9 

1.4257 

1.0 

1.6000 

1.1 

1.7778 

1.2 

1.9592 

1.3 

2.1443 

1.4 

2.3333 

1.5 

2.5263 

X 

X 

1.6 

2.7234 

1.7 

2.9247 

1.8 

3.1304 

1.9 

3.3407 

2.0 

3.5556 

2.1 

3.7753 

2.2 

4.0000 

2.3 

4.2299 

2.4 

4.4651 

2.5 

4.7059 

2.6 

4.9524 

2.7 

5.2048 

2.8 

5.4634 

2.9 

5.7284 

3.0 

6.0000 

\ 


As  in  the  previous  case,  the  coordinate  transformation  given  in  Equation  (118) 
is  used  for  z,  and  no  transformation  is  used  for  y.  The  dimensions  of  the 
physical  domain  in  the  x,  y,  and  z  directions,  respectively,  are  6.0  x  4.0  x 
5.0,  and  the  dimensions  of  the  computational  domain  are  3.0  x  4.0  x  1.5  with 
32,  42,  and  17  nodes  taken  to  establish  the  finite-difference  grid.  The  flow- 
field  was  calculated  for  this  case  with  Re  »  100,  and  the  ground-plane  pres¬ 
sure  distribution  and  x-component-of-velocity  profiles  for  z  -  0  are  shown  in 
Figures  13(a)  and  (b). 

When  the  primary  jets  are  of  unequal  diameter  and  unequal  strength  as 
they  discharge  through  the  blocking  surface,  the  computational  domain  used  to 
calculate  the  flowfield  is  shown  in  Figure  14.  It  is  bounded  by  the  three  jet 
symmetry  planes  (x  *  0,  x  *  S,  and  z  *  0),  the  outflow  plane  (z  =  W),  and  the 
two  solid  surfaces  (y  ■  0  and  y  =  H).  For  each  of  the  jets,  the  entering 
velocity  profile  is  evaluated  using  Equations  (119)  and  (120). 

The  coordinate  transformation  defined  by  Equations  (122)  through  (129)  is 
used  in  the  region  0  <  x  <  S/2  (Scj  =  1.5,  SC2  =  0.5,  Spj  *  2.0,  and  Sp2  = 

1.0)  and  is  reflected  about  the  raidplane  of  the  computational  domain  for  S/2  < 
x  <  S.  The  values  of  x  as  a  function  of  x  are  tabulated  below. 
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Boundary  conditions 


x 

x 

y 


=  0:u  =  0.£=0.£-0, 
-S:«-0.#-0.fS«0. 
=  0:  u  =  0,  v  =  0,  w  =  0,  k 


ik.  =  o  - 

ax  0| 

=  o,  e  * 


Re 


& 

ay2 


y  ■  H:  u  «  0,  v  =  v  (x,  H,  z),  w  »  0,  k  =  k  (x,  H,  z),  e 

z  ■  0*  ^  »  0  ^  s  o  w  *  0  ^  B  q  35  q 
"dz  =  '  dz  '  '  3z  '  8z  “ 

*-  W:lr  =  0--3i*0-lr-0'ii'0'-|f*0 


C  (x,  H,  z) 


GP2l4TTt*12 


Fiptt  14.  Ctmpatatioaal  doaiaia  for  aaeqaatatreagth  jets  with  aomal 
presence  of  a  blocking  plate. 


lathe 


0. 

0. 

2.1 

3.2195 

0.1 

0.1017 

2.2 

3.4286 

0.2 

0.2069 

2.3 

3.6279 

0.3 

0.3158 

2.4 

3.8182 

0.4 

0.4286 

2.5 

4.0000 

0.5 

0.5455 

2.6 

4.1739 

0.6 

0.6667 

2.7 

4.3404 

0.7 

0.7924 

2.8 

4.5000 

0.8 

0.9231 

2.9 

4.6531 

0.9 

1.0588 

3.0 

4.8000 

1.0 

1.2000 

3.1 

4.9412 

1.1 

1.3469 

3.2 

5.0769 

1.2 

1.5000 

3.3 

5.2075 

1.3 

1.6596 

3.4 

5.3333 

1.4 

1.8261 

3.5 

5.4545 

1.5 

2.0000 

3.6 

5.5714 

1.6 

2.1818 

3.7 

5.6842 

1.7 

2.3721 

3.8 

5.7931 

1.8 

2.5714 

3.9 

5.8983 

1.9 

2.7805 

4.0 

6.0000 

2.0 

3.0000 

The  coordinate  transformation  given  in  Equation  (118)  is  used  in  the  z 
direction,  and  no  transformation  is  used  in  the  y  direction.  The  dimensions 
of  the  computational  domain  are  4.0  x  2.0  x  1.5,  with  42,  22,  and  17  grid 
points  used  in  the  x,  y,  and  z  directions,  respectively. 

Two  flowfields  were  calculated  for  the  case  of  two  unequal-diameter  jets 
with  normal  impingement.  In  both  cases  S/Dj  -  6,  H/Dj  -  2,  and  Re  -  100, 
where  the  Reynolds  number  is  based  on  the  exit-plane  properties  of  jet  1.  In 
the  first  configuration  d2/di  ■  0.75  (with  rgj  ■  0.4,  rbj  ■  0.5,  re2  *  0.3, 
and  rb2  -  0.375  in  the  entering  jet  velocity  profiles),  and  in  the  second 
configuration  d2/dj  «  0.515  (rgl  -  0.4,  rbj  -  0.5,  re2  *  0.206,  and  rb2  * 
0.2575).  The  ground-plane  pressure  profiles  and  x-component-of-velocity 
profiles  for  z  -  0  are  shown  for  these  two  cases  in  Figures  15  and  16,  along 
with  the  pressure  data  measured  by  Jenkins  and  Hill,  Reference  20,  for  the 
same  configuration  in  the  absence  of  a  blocking  plate  with  a  Reynolds  number 
on  the  order  of  100  000.  Again,  the  flowfleld  model  provides  the  observed 
pressure  distribution  in  the  impingement  zones  but  fails  to  provide  the 
measured  pressure  variation  in  the  fountain  because  of  the  diffusion  effects 
for  Re  -  100. 
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5.  SUMMARY 


Using  Che  complete  elliptic  conservation  equations  for  incompressible, 
steady,  three-dimensional,  viscous  flow,  a  flowfield  model  has  been  solved 
numerically  for  isolated  impinging  jets  and  for  two  interacting  impinging  jets 
with  fountain  formation.  The  time-averaged  Navier-Stokes  and  turbulence  model 
equations  were  written  in  terms  of  Poisson  equations  for  the  velocity  compo¬ 
nents  and  transport  equations  for  the  vorticity  components,  turbulent  kinetic 
energy,  and  turbulent  dissipation.  Coordinate  transformations  were  introduced 
into  the  system  to  extend  the  physical  solution  domain  in  the  far  field  where 
relatively  simple  constraints  could  be  imposed  on  the  flow  variables  as  bound¬ 
ary  conditions.  The  Poisson  equations  were  discretized  using  the  conventional 
central-difference  algorithm,  and  the  transport  equations  were  discretized 
using  Agarwal's  third-order-accurate  upwind-difference  scheme.  The  finite- 
difference  forms  of  the  Poisson  equations  were  solved  with  the  tridiagonal 
algorithm,  and  the  finite-difference  forms  of  the  transport-type  equations 
were  sol-.ed  with  a  pentadiagonal  matrix-inversion  routine.  Alternating- 
direction  line  iteration  was  used  to  update  the  variables.  Three-dimensional 
flowfields  were  computed  for  laminar  channel  flow  (a  test-case  problem);  iso¬ 
lated,  turbulent,  impinging  jets  with  normal  and  inclined  impingement;  and 
interacting,  turbulent,  impinging  jets  with  fountain  formation. 

Based  on  these  computations,  the  following  conclusions  are  made  with 
regard  to  the  flowfield  model  and  the  solution  scheme  used  in  the  present 
work: 

(1)  Solution  of  the  time-averaged  Navier-Stokes  equations  in  terms  of 
velocity  and  vorticity  is  more  efficient  than  solution  of  the  equations 
in  terms  of  scalar  and  vector  potential  functions  and  vorticity  used  in 
the  previous  jet-impingement  work,  Reference  8. 

(2)  The  Jones-Launder  two-equation  turbulence  model  provides  a  better 
description  of  the  turbulence  field  than  the  Glushko  one-equation  model 
used  in  the  previous  jet-impingement  work.  Reference  8,  since  the  turbu¬ 
lence  length-scale  distribution  need  not  be  specified  throughout  the 
flowfield  but  rather  is  computed  as  part  of  the  solution. 


(3)  The  coordinate  transformations  used  to  specify  the  far-field  boundary 
conditions  provide  a  more  general  solution  scheme  since  they  eliminate 
the  need  to  define  the  boundary  conditions  in  the  near  field  on  the  basis 
of  empirical  data,  as  was  required  in  the  previous  jet-impingement  work, 
Reference  8. 

(4)  Agarwal's  third-order-accurate  upwind-difference  scheme  used  to  discre¬ 
tize  the  transport-type  equations  requires  less  computing  time  and  stor¬ 
age  than  Hoffman's  augmented-central-dif ference  scheme  which  was  used  to 
calculate  the  impinging  jets  described  in  Reference  8. 

(5)  The  use  of  coarse  finite-difference  grids  imposed  by  computer  storage 
constraints  limits  the  Reynolds  number  for  which  convergent  solutions  of 
the  governing  equations  can  be  achieved.  The  surface  pressure  distribu¬ 
tion  computed  for  Re  =  100  compares  well  with  measurements  obtained  at 
higher  Reynolds  numbers  in  the  impingement  region  for  both  isolated  and 
interacting  jets.  However,  in  the  latter  case  the  solutions  for  a 
Reynolds  number  on  the  order  of  100  do  not  provide  the  measured  pressure 
increase  in  the  fountain  for  a  Reynolds  number  on  the  order  of  100  000 
because  of  the  relatively  large  diffusion  effects  in  the  low-Reynolds- 
number  results. 

A  potential  approach  for  extending  the  present  Navier-Stokes  procedure  to 
higher  Reynolds  numbers  within  the  storage  constraints  of  a  CDC  7600  computer 
is  to  employ  wall  functions  in  the  calculations.  In  the  impinging- jet  work 
completed  to  date  at  MDRL,  the  governing  equations  have  been  integrated  to  the 
wall  through  regions  of  steep  velocity  gradients  using  damping  functions  in 
the  turbulence  model  to  provide  the  proper  low-Reynolds-number  behavior  near 
the  solid  surface.  For  Reynolds  numbers  on  the  order  of  10^,  this  approach 
would  require  an  extremely  fine  finite-difference  mesh  near  the  solid  bound¬ 
aries,  with  an  associated  penalty  in  increased  computing  time  as  well  as  in 
increased  storage.  Using  wall  functions,  which  are  described  in  Reference  21 
for  two-dimensional  flows,  the  solutions  of  the  time-averaged  Navier-Stokes 
equations  and  high-Reynolds-nuraber  form  of  the  Jones-Launder  turbulence  model 
equations  are  patched  onto  fully  turbulent,  local  equilibrium  wall-law  pro¬ 
files.  With  this  method,  damping  functions  are  not  required  in  the  turbulence 
model,  and  the  governing  equations  need  not  be  integrated  all  the  way  to  the 
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wall.  The  chief  merit  of  this  approach  is  that  it  economizes  computer  storage 
and  time  because  it  eliminates  the  need  for  a  fine  f inite-dif ference  mesh  near 
the  solid  surface.  The  additional  points  that  would  he  required  near  the 
solid  boundaries  without  the  use  of  wall  functions  can  be  distributed  else¬ 
where  in  the  flowfield. 

Although  application  of  the  time-averaged  Navier-Stok.es  and  turbulence 
model  equations  to  three-dimensional  impinging  jets  is  a  complex  computational 
problem,  the  present  work  has  demonstrated  reasonable  agreement  between  the 
calculated  and  experimental  flow  properties.  The  disparities  between  the  com¬ 
putations  and  measurements  are  not  due  to  a  deficiency  in  the  flowfield  model 
but  rather  to  the  Reynolds-number  limitations  imposed  by  the  storage  capacity 
of  the  CDC  7600  computer.  No  alternate  computational  method  that  does  not 
rely  on  extensive  empirical  data  has  been  reported  to  date  that  can  treat  the 
three-dimensional,  interacting  impinging  jets  as  accurately  as  the  method 
described  in  the  present  work. 
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